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ABSTRACT 


^jVecenl evidence has shown that the NASA/Lewis 
Ice Accretion Model, LEWICE , docs not predict accurate 
ice shapes for certain glaze ice conditions. This paper will 
present the methodology used to make a first attempt at 
improving the ice accretion prediction in these regimes. 
Importance is given to the correlations for heat transfer 
coefficient and ice density, as well as runback flow, selec- 
tion of the transition point, flow field resolution, and drop- 
let trajectory models. Further improvements and 
refinement of these modules will be performed once tests 
in NASA’s Icing Research Tunnel, scheduled for 1993, are 
completed. 


NOMENCLATURE 


A Area (m 2 ) 

A, Area at surface of bead (m 2 ) 

A, Area at top of bead (m 2 ) 

b Bead height (m) 

b, v Rectangular bead height (m) 

C p Heat capacity (J/kg'K) 

e Vapor pressure (N/m 2 ) 

F Wetness factor (dimensionless) 

g Acceleration due to gravity (m/s 2 ) 

h Heat transfer coefficient (W/m 2 'K) 

k Thermal conductivity (W/m*K) 

jC, Lewis number (dimensionless) 

L f Latent heat of fusion (J/kg) 

L, Latent heat of vaporization (J/kg) 

LWC Liquid water content of air (kg/m 3 ) 
m Mass (kg) 

MW Molecular weight (kg/kg mol) 

N f Freezing fraction (dimensionless) 

p Static pressure (N/m 2 ) 

R Radius of sphere (m) 

R 0 Radius of original drop (m) 

R„ Mass flux of impinging water (kg/m 2 s) 
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INTRODUCTION 


r Recovery faclor (dimensionless) 

s Diameter of spread drop (m) 

S Spread faclor (dimensionless) 

T Tempcralurc ('K) 

AT Tempcranirc difference f K) 

l Time (s) 

V Volume (m 5 ) 

V. Right velocity (m/s) 

y Perpendicular to the blade surface (m) 

Greek Letters: 

{} Collection efficiency (dimensionless) 

Viscosity (kg/m-s) 

v Thermal diffusivity (m 2 /s) 

v f Kinematic viscosity (m 2 /s) 

0 Contact angle (radians) 

p Density (kg/m 3 ) 

c Surface tension (N/m) 

x Time scale (sec) 

x f Shear stress (N/m 2 ) 

Subscripts: 

a air 

c convection 

cond conduction 

e edge of boundary layer 

evap evaporation 

f fluid 

kc kinetic energy 

1 liquid 

o total property 

rec recovery property 

s solid 

sens sensible heat transfer 

w water 

oo frce-stream properly 

Abbreviations: 

LEWICE/P Original Two-Dimensional Potential 

Row Ice Accretion Program 

LEW ICE/1 BL Ice Accretion Program Using Potential 

Flow' with Interactive Boundary Layer 
LEWICE/E Euler Row Ice Accretion Program 


LEWICE/NS Navicr-Slokes Row Ice Accretion 

Program 

LEW1CE/T LEWICE/P Program with Conduction 

and Electrothermal Deicer Modelling 
Capability 

LEWICE/TNG LEWICE/P Program with Improved Ice 
Physics Model 

INT RODUCTION 

^^ecendy, several advancements have been made 
in the development of ice accretion theory. It is the pur- 
pose of this report to incorporate these advancements into 
the NASA/LEWIS ice accretion program LEWICE 1 . Spe- 
cifically, attention is focused in two areas. First, a method- 
ology is formulated for predicting the 'sand-grain’ 
roughness model which controls transition and the level of 
convective heat transfer. Second, a new formulation is 
developed for modelling runback water flow. The previous 
methodology, originally formulated by Messinger 2 , results 
in a constant flow of runback water at the surface for glaze 
ice conditions. However, based on the evidence presented 
by Olsen 3 and more recently by Hansman 4 , it has been 
shown that there is a short transient at work in the ice 
accretion process, and that past this initial phase, much 
less surface water flow occurs. This paper will present a 
new runback model in an attempt at modelling this initial 
transient. 

The first section of this paper will be to present 
enhancements in the LEWICE model based on previous 
experimental data and theoretical advancements devel- 
oped by other researchers. It was found that although sev- 
eral researchers had made specific improvements to the 
program, or had performed tests for the purpose of 
improving this program, a single code which combined all 
of these previous achievements did not exist. Specifically, 
improvements made by Cebeci 5 in fluid flow, Rios 6 in ice 
density, Poinsatte 7 in heat transfer, Hansman 4 and 
Yamaguchi 8 in transition modelling, Al-Khalil 9 in runback 
flow and Wright 10 in heat conduction were incorporated 
into this program. 

The second section of this paper will be to develop 
equations which model the transient freezing of a single 
drop of water on the airfoil surface. The equations for the 
geometry of a drop will also be developed. Also, the flow 
rate of a liquid film will be determined using this analysis, 
following the methodology developed by Al-Khalil 9 . As 
droplets freeze, they will form a roughened surface much 
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more rough than the original airfoil surface and much 
rougher than typical ‘sand-grain’ roughness. It is the pur- 
pose of this section to make a first attempt at predicting a 
characteristic roughness height based on a limiting case 
analysis of surface tension controlled water flow. 

IMPROVEMENTS TO THE MODEL 

Flow Field Resolution 

The potential flow solution for iced geometries can be 
very irregular, resulting in several calculated stagnation 
points. Since the integral boundary layer technique used in 
LEWICE needs a single stagnation point to start at, the cal- 
culated iced geometry is smoothed so that only one stag- 
nation point is calculated. The technique used in the 
current version was developed by Cebeci 5 . He calculates a 
smoothed ‘pseudo-surface’ which is used only in the cal- 
culation of the flow field. In the current model, the flow 
solution was smoothed instead of the airfoil coordinates. 

The potential flow code is preferred for the develop- 
ment of a new ice accretion theory over more complex, 
Euler, or Navier-Stokes codes available due to the faster 
execution time. As this ice accretion module become fur- 
ther developed, an analysis will be performed to determine 
the necessary complexity of the flow solution. 

Droplet Trajectories 

The enhancements made to the droplet trajectory 
algorithm are primarily concerned with the speed of the 
calculation. Since approximately 80% of the computa- 
tional time in LEWICE is spent calculating droplet trajec- 
tories, methods of speeding up this module were 
investigated. 

Three improvements were made in this area. First, the 
initial point for calculation was started at a y-Iocalion 
determined by the angle of attack of the flow, not parallel 
to the airfoil. This results in fewer trajectories needed to 
find the desired range. Even at an angle of attack of 4\ the 
program needs to compute only three trajectories to find 
the desired range, whereas previously it required six. 

Second, during the determination of the impingement 
limits, the next starting point for a new trajectory was 
determined using a ‘weighting’ factor. Previously, the next 
starting point for a new trajectory was halfway between a 
missed trajectory and a trajectory which hit the surface. 
Currently, the next trajectory is started between these two 
limits, however, its starting location is ‘weighted’ such 


that the resulting trajectory will travel closer to the top (or 
bottom) of the airfoil. This results in fewer trajectories 
needed to find the impingement limits. This technique has 
not been as successful at decreasing the number of trajec- 
tories. Usually only one or two trajectories out of twelve 
are saved using this technique. 

Finally, methods were investigated to reduce the num- 
ber of air velocities needed. During each time step of each 
trajectory, the algorithm performs a predictor-corrector 
iteration to solve the integration of the non-linear momen- 
tum equation. For this integration, the program requires 
the air velocity at that location. Previously, this has been 
found by summing the contribution from each panel in the 
potential flow solution. So far, two methods have been 
attempted to simplify this cumbersome computation. 

First, the air velocity was simply lagged one time 
step. Thus for each iteration in the predictor-corrector, the 
same value of the air velocity was used. This assumption 
is valid, as the difference in position between predictor 
and corrector is much smaller than the distance between 
the drop position at time step n and its position at time step 
n + 1. This modification in itself can decrease the compu- 
tation time by a factor of three and thus merits further 
investigation in three-dimensional codes as well as more 
complex two-dimensional codes. 

Second, instead of calculating the air velocity for each 
step in the trajectory routine, velocities were calculated on 
an off-body rectangular grid and an interpolation scheme 
was used to find values in between grid points. This results 
in fewer air velocities which need to be calculated. This 
method also reduced the computation time, although the 
improvement is not as great if a fine mesh is used. This 
happens because the number of panels is small in a two- 
dimensional code. Further investigation is warranted in 
three dimensional codes for this method. 

Finally, one correction was made to the droplet trajec- 
tory algorithm to increase its accuracy. This correction is 
concerned with the ‘d-shiff used in LEWICE . When cal- 
culating the off-body air velocity at any location using a 
potential flow solution, the answer will not be reliable if 
the point lies too close to any given panel. To alleviate this 
problem, LEWICE creates a set of coordinates which lie a 
distance *d-shift* normal from the surface. 

Previously, when a drop hits the enlarged surface, it 
was treated as if it hit the actual surface. This results in an 
over prediction of the collection efficiency because drops 
which might impinge the larger d-shifted* surface might 
not hit the actual surface. Currently, the drops arc allowed 
to proceed until they impinge the actual surface. However, 
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if ihc distance between the drop and a given panel is less 
than *d-shift\ the air velocity which is used will be inter- 
polated from the value at a location ‘d-shift* from the 
panel and the value at the panel which is provided from 
the potential flow solution. This is illustrated in Figure L 



Figure 1. Location of Drop Within ‘D-Shift 1 of Surface 


Heat Transfer Coefficient 

The heat transfer coefficient in LEWICE is controlled 
via the skin-friction which in turn is controlled by the 
equivalent sand-grain roughness. In glaze ice conditions, 
when all of the ice does not freeze upon impact, the local 
freezing rate can be greatly controlled by the heat transfer 
coefficient The current model still relys on the sand-grain 
roughness; however, the correlation itself has been cor- 
rected to more closely reflect heat transfer coefficient data 
taken by Poinsatte 7 using the NASA/Lewis Twin Otter 
Research Aircraft and using the Icing Research Tunnel. 

In that study, hemispherical elements were glued onto 
a NACA0012 airfoil with heater elements in the airfoil. 
The heaters were kept at a constant temperature of 60 °F 
and the necessary heat supplied was recorded for various 
Reynolds numbers and angles of attack. Four different 
roughness patterns, including a clean airfoil, were used in 
this study. 

These tests are not considered ideal for use in 
LEWICE , as it does not use actual ice roughness. How- 
ever, they provide a means of modifying the correlation 
for large sand grain roughnesses. Since the element 
heights were 1 mm, a sand grain roughness of 1 mm was 
used for comparison with the LEWICE correlation. 
LEWICE greatly overpredictcd the heal transfer coefficient 
using this value. Three corrections were made which 
resulted in better agreement with the experimental results. 


First, the temperature dependance of both the thermal 
conductivity and viscosity were corrected. Based on what 
appears to be a typographical mistake in the coding, the 
temperature dependance of these variables was too high. 

Second, the computation of the transition Reynolds 
number was changed when the roughness height exceeds 
the boundary layer thickness. Transition will still occur 
when the Reynolds number based on roughness exceeds 
600. However, if the roughness height exceeds the bound- 
ary layer thickness, the definition of the Reynolds number 
is changed. Previously, the velocity was set to the velocity 
at distance 5, i.e., the potential flow solution. The rough- 
ness height was still used in the calculation. Currently, the 
routine uses the boundary layer thickness instead of the 
roughness height in the Reynolds number calculation. 


Previous Definition: Re k 


v* K 

v 


for kj > 5 


( 1 ) 


F 6 5 

Current Definition: Re k = for k, > 5 


( 2 ) 


Third, the turbulent roughness value was calculated 
by the roughness height model described later in this 
paper. Any roughness in excess of the turbulent boundary 
layer thickness is treated as follows: 

1) The roughness height up to the boundary layer 
thickness direedy affects the convective heat transfer as a 
sand-grain roughness. 

2) The amount beyond this level affects the convec- 
tive heat transfer only due to the flow acceleration over 
this protuberance which can be calculated by the flow 
solver. 

Using these corrections, much better agreement was 
achieved with Poinsatte’s experimental data. Future tests 
are planned which will attempt to measure heat transfer 
coefficient on plastic ice models. 


Ice Density 

The ice density correlation used in LEWICE was 
developed by Macklin 11 using tests on rotating cylinders. 
The Macklin correlation is based on three parameters: the 
local surface temperature, the volume median droplet 
diameter, and the droplet impact velocity. This correlation 
is believed to be inaccurate because it is based on a param- 
eter which was not measured in the test (droplet impact 
velocity). Macklin calculated the impact velocity from 
measured variables using a very simplified analysis. 
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Rios'* presented an alternative ice density correlation 
hased on tests performed by Jones 12 . These tests were also 
performed on rotating cylinders and hence arc not ideal for 
ice prediction on airfoils. It is based on five parameters; air 
temperature, air velocity, liquid water content, volumetric 
mean droplet diameter, and cylinder diameter. This rela- 
tionship can, at best, only provide the average ice density 
as it is based on ambient parameters. 

The implementation in LEWICE assumes that the 
functional form of the correlation can be applied at each 
chord location using the local equivalent to the ambient 
variable. The surface temperature, collection efficiency, 
and local radius of curvature are substituted into the corre- 
lation. This results in a chordwise variation of ice density 
for most conditions, with the higher glaze ice density at 
the leading edge, progressing to lower rime values down- 
stream. The form of this correlation is given below. 

/«( p lW ) = -0.15 x (1 +6043.T 265 ) (3) 

where 

d° n *]x V^ 59 x (PLV^C) 021 

S O (~T air ) 0 2i ~ (4) 

Future tests are planned which will measure the local den- 
sity on airfoil ice shapes to verify this correlation. 

Runback Water Flow 

The methodology used in LEWICE assumes that 
water which impinges and does not freeze will automati- 
cally flow into the next control volume. This can result in 
relatively large water flow rates on the surface for glaze 
icc conditions. These large flow rates however arc not sup- 
ported by qualitative assessment of icc accretion tests or 
by close up videos. 

The current model uses a formulation derived by Al- 
Khalil 9 for runback water flow. In his model, water flow is 
assumed to be shear driven by the air flow. His formula- 
tion solved for the complex, two-dimensional flow of sur- 
face water, typically using engine inlets equipped with an 
anti-icer. Since LEWICE is a two-dimensional flow pro- 
gram the airfoil surface is one-dimensional, hence the the- 
ory developed by Al- Khali 19 can be simplified in this 
paper. 

This model has been further modified by assuming 
that the ‘wetness factor* described by Al-Khalil can be 
independently determined by the amount a drop spreads 


upon impact. This is determined from experimental data 
by Macklin H and by Hansman 4 who measured the contact 
angle of the drop from which the spread factor can be 
determined. Additional tests are needed to confirm the 
relationship between contact angle and the ambient param- 
eters. 

Additionally, the water is allowed to stagnate if the 
force of the flow is found to be less than the surface ten- 
sion force. This model has the advantage of modelling 
standing water in each control volume, and should be 
accurate for modelling separation regions in Navidr- 
Stokes flow, since water will not flow into a separation 
region. 

Heat Conduction 

Wright 10 modelled the two-dimensional transient con- 
duction in an airfoil with ice accretion for the application 
of modelling an electrothermal deicer. When all internal 
heat sources are turned off, this program will predict the 
ice accretion on an airfoil with conduction effects. How- 
ever, this simplified application can be handled without the 
additional computational burden. An analytic solution is 
available which gives a reasonable approximation for this 
case. Locally, at each chordwise location, the heat conduc- 
tion can be modelled using a semi-infinite flat plate 
assumption where the airfoil surface al time = 0 is sud- 
denly raised from the recovery temperature to the icing 
temperature (normally 0*C in glaze ice conditions). This 
assumption has been verified by performing parameter 
runs using the full conduction model. Since the airfoil skin 
thickness is not infinite, this assumption can remain in 
effect only until the heat has penetrated through the airfoil 
skin. The parameter used to measure this is termed the 
penetration thickness. When the penetration thickness 
equals the airfoil thickness, the heat has penetrated 
through the airfoil skin. At this point in time, the conduc- 
tive flux at the surface is small enough to be ignored. 

Liquid Water Content 

The physical evidence obtained from tests performed 
so far show that the final ice shape is greatly dependant 
upon the physics during the first few seconds. It is felt that 
a large amount of the repeatability problems in generating 
experimental ice shapes stems from an inability to exactly 
match, either in flight or in the IRT, the exact variation in 
conditions. This is especially true during this initial time 
frame. It has been established in the IRT that liquid water 
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content will increase from 0 to the desired LWC value 
within about 20 sec. The program currently is capable of 
modelling this transient linearly. At present, no parametric 
studies have been performed to determine the effect of dif- 
ferent LWC transients on the final ice shape, nor have any 
flight cases been evaluated for this purpose. 

GLAZE ICE ACCRETION MODEL 


As stated earlier, the experimental evidence as shown 
by Olsen 3 and Hansman 4 clearly show that the current 
methodology for ice accretion is incorrect when applied to 
glaze ice. However, much of this data is qualitative in 
nature, which makes it difficult to incorporate into a com- 
puter model. The following analysis presents a first step at 
modelling some of these phenomena. Where quantitative 
data exists, this model will be compared to these experi- 
ments. 

Bead Geometry 

Assume that as a spherical droplet of water impinges 
a locally flat surface and deforms, it retains a spherical 
form with a section cut from the bottom (see Figure 2). 



Figure 2. Dispersement of Water Droplet 


Furthermore, assume that in the short time it takes this 
drop to stabilize that the volume remains constant. This 
assumes that evaporation losses are small and that little of 
the ice has frozen. As it is assumed the drop takes on the 
order of 10‘ 7 sec. to stabilize at any given contact angle, 
these assumptions are reasonable. This lime scale is 
derived by dividing the drop diameter by the drop velocity. 
This volume at any given angle can be determined as the 
volume of a sphere with a conical section cut from it plus 


the volume of a fiat-bottomed 
This is represented in Figure 3. 



cone with an equal angle. 



r 


Figure 3. Representation of Partial Volume of Sphere 


This volume can be represented mathematically by 
adding the volume of a sphere with a conical section cut 
out plus the volume of a flat-bottomed circular cross-sec- 
tion cone. The first of these volumes is found by integrat- 
ing 


R * K 

V = JJJrcr 2 sin0sin<()^<{>d0dr (5) 

0 0 a 

if this equation were integrated from 0 to k instead of a to 
7r, the equation for the volume of a sphere is obtained. 

The volume of a cone is given by 

it i 

v = yh ( 6 ) 

or by writing this in terms of the sphere diameter R 
(the hypotenuse of the triangle with sides h and r), 

V = ^/? 3 (cosa-(cosa) 3 ) (7) 

The total volume of this section is then 

V = ^# 3 (2 + 3cosa-(cosa) 3 ) (8) 

where again R is the radius of the sphere, and a is the 

angle at which the sphere is cut off (see Figure 3 above). 

For a > it! 2, the segment of the sphere which contains 
the mass of the drop is less than half of the total volume of 
that sphere. The volume can be represented by subtracting 
the above volume from the total volume. This yields 

V = (2-3cos0 + (cos0) 3 ) (9) 
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However, this equation gives the volume in terms of 
an angle which is equal to n - a. To be consistent, the 
angle can then be changed back to the same point as the 
above volume (see Figure 4). 



Performing this substitution yields 

V = ^/? 3 (2 + 3cosa- (cosa) 3 ) (10) 

which is the same equation obtained before for a < rc/2, 
hence the above equation describes the deforming sphere 
at any angle. It is assumed that as more drops impinge or 
as drops coalesce that only the mass is changed and that 
the geometry remains the same. 

Energy Equation 




L r M w h r T t 

+ RS P ,(T„ - TJ + (<*, - e s y) 

'o, * 


- R * L / N / (ii) 

T rec represents the recovery temperature, as defined by 
Schlichting 14 . 

However, this analysis does not account for the differ- 
ence in area between the icing surface (which is presumed 
flat) and the surface area of the bead. Furthermore, it is 
more convenient for the current analysis to replace the 
accretion rale R w N f as p.db/dt, where instead of using the 
freezing fraction, the increasing thickness of the ice db/dt 
is used instead. If A, /A s represents the ratio of the top sur- 
face area to the flat surface area, the energy equation 
becomes 


dT A, i 


L v M Jc A , T e 

+ *„C pi (T„-T m ) + -L-l-E-i (,, - , -') 


P » C P . A s 


The energy equation for a given drop at known con- 
tact angle is found by making the following assumptions: 

' 1) No part of the drop goes below freezing until the 

entire drop is frozen, hence = T TOlt ; 

2) Conduction heat loss and sensible heal transfer 
occurs only at the interface between the surface and the 
drop; 

3) Evaporation and convection occurs at the water/air 
interface; 

4) Ice forms at the bottom surface of the drop; and. 

5) The thermal properties of the solid surface can be 
found by adding the thermal resistances of each layer 
(including the ice) in a composite structure. 

Following the derivation of Wright 10 , the energy bal- 
ance on an individual drop is given by 



( 12 ) 


Since the drop is very small compared to the thickness 
of the surface, and since the drop temperature can be 
assumed to be close to T^,, until it freezes completely, the 
heat loss into the surface can be expressed as the analytic 
solution of heating a semi-infinite slab, hence 


k dr\ _ k ( T rec- T J 
dy y = o Jnvr 


(13) 


where v is used as the thermal diffusivity in order to avoid 
confusion with the angle a. To simplify future equations, 
the energy equation is rewritten in terms of a temperature 
difference times an effective h for that term. This yields 
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T 05 A 

^ cond ( AT rec) ( 7 ) + 


where 


db 


+ Ktns ( A + A T fv ) = 


(14) 


l cond 


Jkvx 


h sens ~ 


A Tt ‘ 2 C + T ' ec T ~ 

Pt 


AT = T ~T 

rec m * rec 


A7\ 




evap 


r a M a c p /? 




(15) 

(16) 

(17) 

(18) 
(19) 


and where t = representative time scale. 

The maximum value of the rate of freezing will occur 
when the rate of freezing above exceeds the rate of 
impingement, which is 


db _ 
dt ~ p s 


( 20 ) 


The amount of water on the surface at any time can be 
found by subtracting the rate of ice formation from the rate 
of impingement and integrating. Note that since the quan- 
tity R w is defined in terms of the liquid water content, this 
value is not necessarily constant. 

The following geometric parameters can now be 
defined as follows. Let b iV be defined the average height of 
the bead, s is the diameter of the spread bead, and b is the 
bead height, as shown in Figure 5. 



Figure 5. Geometry of Water/Ice Bead 


= 2/? sin 9 

(21) 

R( 1 - cos0) 

(22) 


where 0 is the contact angle, which is defined as the tan- 
gent angle of the drop where it intersects the surface. Note 
that 0 = 71 - a. The volume of the bead can be expressed as 
the volume of a cylinder with diameter s and height b. v or 
as the partial volume of a sphere as before. Equaling these 
two yields 

£ s 2 b av = ti/? 3 (2- cos0 + (cos0) 3 ) (23) 

Substituting the relations for b and s into this expression 
yields 


(2+cos0) 
flV (1 + COS0) 


(24) 


The area ratio is the ratio of the exposed surface area 
to the area at the base. 


A, _ 2nR 2 (l - cos0) 

A s kR 2 ( 1 - cos0) ( 1 + cos0) 


= 2 
(1 + COS0) 

hence the total ice accretion rate becomes 


(25) 


h c (&T rcc + AT) 


evap ' (1 + COS0) 


+ T rec ) 


T 0.5 

( 7 ) +*s,ns(*T, <c + AT„) 



(26) 


Water Flow Criteria 

At a certain point, if enough water forms on the sur- 
face, it will form a continuous film and start to flow. This 
will happen when the shear force on the water overcomes 
the surface tension forces. From Al- Khalii\ 


From simple trigonometry relations. 


dm / _ Fx / h : 
di IVj 


(27) 
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where m f is the mass per unit span, x ( is the shear stress, v 
is the kinematic viscosity, and F is the wetness factor, the 
fraction of the surface covered with water. The velocity of 
the water flow is given by 


Fx / 

v ' - 4" 


(28) 


The force per unit span is obtained by multiplying the 
velocity by the mass flow rate to get 


This quantity has the same units as surface tension, o. 
When this quantity is greater than the surface tension, the 
bead will flow. Performing this operation yields the flow 
criteria for water runback 


1 



where a = surface tension. 

Additionally, the flow can be affected by gravitational 
forces. The method in which gravity is modeled is to 
assume that there is a limiting drop height above which the 
drop becomes unstable at that contact angle and starts to 
flow. If the drops are small, they retain their shape, but if 
more water is added, the drop will reach a mass beyond 
which flow will occur. This will happen when 

mgbcost|r>a4 (31) 

where \|/ is the local angle the airfoil makes with a 
horizontal line. By substituting the proper geometric 
equivalents. 


b> 


110(1 - COS0) 
pgcosxj/ 


(32) 


One of the unknown parameters in the above develop- 
ment is the wetness factor, F (see equation 30). The fol- 
lowing strategy will be adopted here for calculating this 
variable. There is a quantity called the ‘spread factor’, 
which was defined by Macklin 13 as the ratio of the diame- 
ter of the spread drop to the initial drop diameter, and is 
determined by equating the expression for the volume of a 
spread drop with the original spherical volume. 


*/? 3 (2+cosG) (l -cos0) J = *nR y 0 

Substituting the previous expression for the 
drop diameter s into this equation yields 

l 

_ s _ r 4sin0(l -cos8) ~I3 
2 R 0 L(2+ cos0) (1 - cos0)^ 

According to Macklin 13 , a surface is considered com- 
pletely wetted when the contact angle is below 10\ The 
wetness factor can then be defined as the ratio of the 
spread factor at the contact angle 0 to the spread factor at 
an angle of 10\ with an upper limit of F = 1. 

As a check on the previous equation, both Macklin 13 
and Hansman 4 performed experiments to determine the 
variation of drop geometry with temperature. Hansman 
measured the contact angle from close-up movies while 
Macklin measured the diameter of individual frozen drops 
with a micrometer, then divided by the known diameter of 
the incoming drops to obtain the ‘spread factor’. A com- 
parison of the data by these two researchers is presented 
below. 


<331 

spread 


(34) 


! 

bk 



Figure 6. Comparison of Experimental ‘Spread Factor’ 


The upper line in the above plot gives the largest 
spread factor for that temperature while the lower line 
gives the smallest value for that temperature. The differ- 
ence between the upper and lower values represents the 
scatter in the data. The theoretical lower limit shown in 
this plot is computed by inputting the value for a hemi- 
spherical bead which is 0 = n/2. The agreement between 
the two tests show the same trends, however both parame- 
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RESULTS AND DISCUSSION 


icrs arc difficull lo measure accurately, as shown by the 
data scatter. 

Computational Strategy 

The bead height model used in this paper only identi- 
fies the maximum possible roughness level that can be 
obtained. These maximum values are based on physical 
criteria. The actual roughness of the surface is unknown, 
as it could be less than these values. For this analysis, the 
roughness heights used was the smallest of these maxi- 
mums. These criteria are as follows: 

1) Roughness height can not exceed the ice thickness; 

2) Water height can not exceed the height needed to 
flow based on the shear forces exceeding the 
surface tension forces; 

3) Water height can not exceed the height needed to 
flow based on the gravitational forces exceeding 
the surface tension forces; 

4) For both 2) and 3) the total bead height is 
determined using the water height and the freezing 
fraction; 

5) Heal transfer coefficients arc affected by bead 
height only up to the boundary layer thickness, 
whether this is for transition or for turbulent 
values. The additional effect of heat transfer due to 
these beads is presumed to be handled by increased 
velocity in the flow solution. 

The computational strategy for glaze ice accretion is 
as follows: 

1) using a small time step 1 , water is allowed to 
impinge upon the airfoil; 

2) the amount frozen, the liquid bead height and 
runback water flow (if any) are determined from 
the above equations; 

3) the roughness height is set equal to the bead 
height; 

4) boundary layer quantities are recomputed to 
determine new values for heat transfer 
coefficient; 

5) time is incremented; 

6) if the flow is undisturbed 2 by the ice shape, return 
to step 1); 

7) if the flow is changed, compute new flow field 


I The size of this step needed to ensure accuracy has not yet 

been established 

2. The criteria for this has not been established. 


and trajectories, then return to step 1); 

8) if the total simulation time has been reached, stop 
the simulation. 

RESULTS AND DISCUSSION 


The model presented in this paper has an implied 
assumption that the ice accretion process as a whole will 
be modeled better when each of the sub-models presented 
is accurate. Therefore, the first step in verifying this model 
will be to verify, if possible, the sub-model. However, for 
many of these sub-models, this can not be done as yet. 
Therefore, only qualitative assessments can be made as to 
their improvement. Except where otherwise noted, the 
cases used to present the capabilities of the model are 
taken from experimental runs by Shin 15 . Most plots were 
obtained using the following conditions: V„, = 230 mph; 
T 0 = 28 'F; LWC = 0.55 g/m 3 ; MVD = 20 fim; t = 7 min; a 
= 4* 

Flow Field Resolution 

Figure 7 shows a comparison of the output from the 
potential flow solver with the smoothed values after six 
minutes of ice accretion. Even though this solution is more 
uniform, it does not serve its full purpose of determining 
the stagnation point for the boundary layer routine. The 
stagnation point is determined by picking the point or 
points where the velocity changes sign, which still hap- 
pens more than once. At this point, the program selects the 
new stagnation point as the one which is closest to the old 
stagnation point. 

Droplet Trajectories 

There were several modifications made to the droplet 
trajectory routine. Figure 8 shows the optimization of the 
range routine. For this particular case, placing the first x,y 
point at an angle a below the airfoil results in three fewer 
trajectories. 

Figure 9 shows the effect of using a grid-based 
scheme to interpolate velocities. As can be seen in this 
plot, there is an insignificant difference in the path of the 
trajectory using this routine. However, the computational 
savings in 2D were also minimal. The use of this approxi- 
mation will be more beneficial for 3D algorithms. 

Figure 10 illustrates the approximation made when 
off-body air velocities are calculated only at the start of the 
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time step. This plot shows the x,y location of the drop at 
lime n and at lime n+1 every time the velocity routine is 
called. Although it does not show well on this plot, the 
velocity routine is called six times at time n and eight 
limes at time n+1. Currently, it is called once for each 
location. Again, there is no change in the trajectory of the 
drop. As the locations are so close, their individual loca- 
tions can not be seen on the plot. 

Figure 1 1 shows the effect on collection efficiency 
using the interpolation procedure within ‘d-shiFt’ of the 
surface. These results also incorporate the approximations 
described above. As can be seen from this plot, the effect 
of using different ‘d-shifts’ has been greatly minimized. 
Figure 12 shows the difference in the 7 minute ice shape 
and the experimental comparison. Again, very minimal 
difference in ice shape prediction is seen. By contrast. Fig- 
ure 13 shows the variation in predicted ice shape using the 
same two ‘d-shifts*. 

Figure 14 shows the variation in collection efficiency 
with time. Qualitatively, this result appears correct. The 
horns collect more water than the same location on the 
clean airfoil, while fewer drops collect at the stagnation 
point. However, no experimental validation of the quanti- 
tative values have been made, except for a clean airfoil. 
This could be a large source of the error in the predicted 
results. 

Heat Transfer Coefficient 

Figure 15 shows a sample comparison of both the cur- 
rent correlation and the previous correlation with the 
experimental data. As can be seen in this plot, much better 
agreement is achieved using the corrections described ear- 
lier. Since this experiment had only metal roughness ele- 
ments and not ice, the bead height is known. This height is 
used unless it exceeds the boundary layer in either the 
laminar or turbulent regimes. This has the effect of moving 
transition downstream. Assuming that roughness levels 
exceeding the turbulent boundary layer can not be mea- 
sured by the current correlation produced heat transfer lev- 
els approaching that of the experimental data. 

Heat Conduction 

Figure 16 shows a comparison between the approxi- 
mate equation used for heat losses by conduction and the 
computational heat losses by conduction using LEWtCEfr \ 
Good agreement is shown, validating that this approxima- 
tion is reasonably accurate for including conduction 


effects in ice accretion. The variation in the computational 
results arc due to the fact that the heat conduction program 
is set up to analyze electrothermal heaters which will pro- 
duce heat fluxes two orders of magnitude greater than 
these levels. Hence the truncation of the finite-difference 
scheme comes into play for these very low heat flux levels. 

Bead Height Model 

Figure 17 shows the predicted maximum bead heights 
which are used as roughness heights in the program for 
two different time steps. The initial values vary gradually 
between 0.1 mm and 0.8 mm. while at t=360 sec. the vari- 
ation is much more irregular. This is due to the irregularity 
in the potential flow solution. The flow field now has con- 
trol over not only the bead height, but the runback flow as 
well. This advancement will hopefully allow belter flow 
solvers, such as Navi£r-Stoke$ to improve their ice accre- 
tion prediction much more so than a potential flow solver 
because the potential flow solver is not accurately predict- 
ing the flow behavior for glaze ice shapes. 

Ice Density 

Figure 18 shows an example of the chord wise density 
which can be achieved using the current correlation. This 
case was run using a total temperature of -1ST instead of 
the 28*F temperature used for the other cases. This change 
was performed because ice density did not vary for the 
very warm glaze condition ran earlier. This is supported by 
both Macklin’s and Jones’s data which give an average ice 
density approaching the glaze ice value as temperature 
increases. Even though this is a pure rime case, much of 
the ice which forms has a density within 10% of the glaze 
value. This is why rime ice shapes have in the past com- 
pared well with experiments even though the ice shape 
formed in the experiment may have had a lower density. 

Ice Accretion Comparison 

Figure 19 shows a comparison between the predicted 
ice shape and the experimental ice shape presented by 
Shin 15 . He compared experimental ice shapes with those 
predicted by LEW1CE/1BL. Most of these predictions are 
very good. The particular case chosen for comparison here 
was selected because the ice shape predicted by Shin 15 
using LEWICEUBL did not match the experimental data. 
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CONCLUSIONS 


The program created for this paper has been called 
LEWICE/TNG ( LEWICE : The Next Generation) for com- 
parison with previous models. Currently, the overall ice 
prediction with LEWICEITNG is not very much improved 
over LEWICE! P as shown in Fig. 20. It should be noted 
that for this comparison, a roughness height of .5 mm and 
a *d-shift’ of 2% of the chord was used in LEWICE!?. 

The main improvement of the current model is in the 
reduced effect of these immeasurable parameters on the 
solution. Figure 13 showed the variation of ice shape pre- 
diction by LEWICE!? using *d-shifts’ of 2% and 4% of 
chord, while Fig. 2 1 shows the effect of varying roughness 
for the values of 0.1 mm, 0.5 mm, and 1.0 mm using 
LEWICE!?. By contrast, the same exercise was performed 
using LEWICE!TNG. These results are shown in Figs. 21 
and 22. As can be seen in these plots, much less variation 
in ice shape prediction is exhibited using LEWICE/TNG. 
The current variation due to ‘d-shift’ is due to the fact that 
the velocity is not changing linearly near the surface, 
hence using different ‘d-shifts’ in the interpolation proce- 
dure causes a slight difference in the collection efficiency, 
especially near the impingement limits. Since LEWICE/ 
TNG overrides the roughness level input and uses the cri- 
teria described earlier, the input roughness has almost no 
effect on the solution. The difference shown is because if 
there is ice forming at a specific location at time step n+I 
when no ice or water was present at time step n, then the 
model does not have a roughness value available. Until 
this loophole is closed, the input roughness is used in this 
region, hence the slight difference at the icing limits. 


Future Work 

The model used for predicting a roughness height is in 
the early stages of development and needs to be better 
defined. The main benefit of this model is that it allows the 
flow solution to exert greater control over the ice accretion 
prediction. This may not be a beneficial feature when 
using potential flow, as the solution for iced airfoils is very 
irregular. An attempt will be made to incorporate this new 
algorithm into the LEWICE/NS program. The geometry 
modification routines will also be examined. Recently, 
Hansman 4 has suggested alternative ice growth mecha- 
nisms which need to be explored. 

After closely examining the close-up videos, one 
physical effect stands out which is not well predicted by 
this program. This is in the area of collection efficiency 


prediction. Experimental work in quantifying collection 
efficiency has always been performed on dean airfoils. 
Since potential flow programs arc reasonably reliable for 
this situation, the collection efficiency should be well pre- 
dicted. In glaze conditions where horns are produced and 
leading edge separation may occur, the collection effi- 
ciency may not be well predicted. No experimental data is 
yet available to confirm this, so an analysis will be per- 
formed using Navidr-Stokes and potential flow programs 
to identify the difference in prediction of collection effi- 
ciency. The ice shape prediction shown in Figure 19 has 
approximately 15% more ice than the experimental ice 
shape tracing. However, in the experiment, all of the water 
which impinges will freeze somewhere on the airfoil 
(except for a very small amount of evaporation). In the 
program, there is a large amount of runback water flow off 
the airfoil. If all of this runback water were to freeze, there 
would be approximately twice as much ice on the airfoil. 
This means that impingement levels for horn ice shapes is 
well overpredicted by a potential flow program. 

CONCLUSIONS 


Several additional physics of ice accretion were added 
to the existing model. Slight improvement has been made 
to the ice prediction capabilities as a result. The major 
improvement has been a marked decrease in the sensitivity 
of the program to the input values for ‘sand-grain’ rough- 
ness and for ‘d-shift \ Optimization of the droplet trajec- 
tory routines results in a 30%decrease in the overall 
computational time. 
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Figure 7. Effect of Smoothing the Flow Solution 



Figure 8. Optimization of Range Routine 
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Figure 9. Drop Trajectory Using Interpolated Velocities Figure 11. Variation of Bela with D-Shift: LEWICE/TNG 




Figure 10. Drop Position at Two Different Time Steps 


Figure 12. Variation of Ice Shape with D-Shift: LEWICE/ 
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Figure 13. Variation of Ice Shape with d-Shift: LEWICE/P Figure 15. Comparison of FrOssling Number with 

Poinsatte Data 



Figure 14. Change in Collection Efficiency with Time 


Figure 16. Comparison of LEWICE/TNG Conduction 
Approximation to LEWICE/T 
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Figure 17. Chang e in Bead Height Prediction with Time Figure 19. Comparison of LEWICE/TNG Ice Shape with 

Experimental Data 



Figure 18. Sample Variation of Ice Density with Distance 

Figure 20. Ice Shape Comparison for LEWICE/P and 
LEWICE/TNG 
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Figure 21 . Variation of Ice Shape with Roughness: 
LEWICE/P 



Figure 22. Variation of Ice Shape with Roughness: 
LEWICE/TNG 
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